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ABSTRACT 

Goldreich-Sridhar model of incompressible turbulence provides an elegant approach to describing 
strong MHD turbulence. It relies on the fact that interacting Alfvenic waves are independent and 
have random polarization. However, in case of strong interaction a spontaneous local assymetry can 
arise. We used direct numerical simulations to show that polarization alignment occurs and it grows 
larger at smaller scales. This effect would lead to a shallower spectrum and stronger anisotropy. Even 
small changes in these two properties will have important astrophysical consequencies, e.g. for the 
cosmic ray physics. 

Subject headings: turbulence: incompressible, turbulence: magnetic 


1. INTRODUCTION 

Astrophysical plasmas is turbulent and magnetized. 
Magnetohydrodynamic turbulence affects many phenom¬ 
ena including the propagation and acceleration of cos¬ 
mic rays (see Shlickeiser 2003). While foundations of 
the theory of anisotropic MHD turbulence model can be 
traced back to works in 80s (see Shebalin, Matthaeus 
& Montgomery 1983, Higdon 1984, Matthaeus & Brown 
1988) a substantial recent progress was initiated by the 
pioneering study by Goldreich & Sridhar (1995, hence¬ 
forth GS95). There a concept of balancing of linear 
and non-linear term in incompressible MHD equations 
was suggested, which results in a prediction between the 
wavenumbers parallel || and perpendicular T to the local 

direction of magnetic field fey ^ kjj . GS95 model also 

predicts Kolmogorov E{k) ^ spectrum, which is 

consistent with both interstellar (see Armstrong, Rick- 
ett & Spangler 1995) and Solar wind (see Horbury 1999) 
data. 

While the core concept of the GS95 model, namely, 
critical balance, was confirmed in numerical simulations 
(Cho & Vishniac 2000, Maron & Goldreich 2001, Gho, 
Lazarian & Vishniac 2002, henceforth CLV02), the sim¬ 
ulations revealed a difference in the spectral indexes ob¬ 
tained. For instance, Maron & Goldreich (2001) obtained 
the spectral index that is close to —3/2 in contrast to 
—5/3 in CLV02. Index of —3/2 for field-perpendicular 
energy was obtained in simulations of Muller & Grappin 
(2005). This stimulated theoretical efforts to understand 
the difference between the theory and simulations on one 
hand and between different sets of simulations on the 
other hand. In particular, Boldyrev (2005) proposed a 
model in which the interactions are weakened compared 
to the GS95 predictions as a result of 3D structure of ed¬ 
dies. An alternative point of view, namely, flattening of 
the spectrum as the result of intermittency was discussed 
in Maron & Goldreich (2001). 

We note, that the problem of MHD turbulence spec¬ 
trum is of great practical importance. For instance, the 
introduction of GS95 spectrum and its extension to com¬ 
pressible media (see Lithwick & Goldreich 2001, Cho & 
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Lazarian 2002, 2003) resulted in a substnatial changes 
in understanding of cosmic ray propagation and acceler¬ 
ation (see Chandran 2000, Yan & Lazarian 2002, 2004, 
Cho & Lazarian 2006). 

In what follows we address the issue of interaction 
efficiency using numerical simulations. We use both 
compressible and incompressible MHD runs. Contrary 
to some earlier claims, Cho & Lazarian (2002) demon¬ 
strated weak coupling of the compressible and Alfvenic 
parts of the spectra, which justifies comparing Alfvenic 
turbulence in incompressible and compressible media. 

2. INCOMPRESSIBLE MHD 

In the incompressible case it is convenient to use 
Elsasser variables w = v^-|-v—b, u = —v^-|-v-|-b, where 
Vyi and b are mean and total magnetic fields in velocity 
units (see Biskamp 2003). With those MHD equations 
will have a symmetric form 

9tU — (v^ • V)u = —(w • V)u — VP, (1) 

dt'w + (vA • V)w = —(u • V)w — VP, (2) 

V • u = V • w = 0. (3) 

Here the total pressure P = p + is determined 

by incompressibility conditions (3). Explicitly 

d^x' Vw : Vu 
471 |x' — x| 

Aside from symmetry, the remarkable property of these 
equations is the existence of the exact nonlinear solutions 
for one field, such as u = f(r -|- v^t), in the absence of 
the other, w = 0. 

Variables u and w can be Fourier decomposed into 
waves which have a dispersion relation of a; = ra^|| • 
These are actually consist of two modes, shear Alfven 
waves, which have u and w perpendicular to both k 
and va, and pseudo-Alfven, or slow waves that in in¬ 
compressible fluid compress magnetic field. It have been 
shown both euristically and numerically, that in devel¬ 
oped turbulence shear waves govern the cascade and the 
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back-reaction of slow waves are relatively unimportant 
(see Maron & Goldreich 2001). 

In equations (1) and (2), written in Fourier space, lin¬ 
ear term could be estimated as VAk\\u and nonlinear as 
wk±u. With sufficiently small wave amplitudes linear 
term dominates and the theory of so called weak or wave 
turbulence could be built (Galtier et al, 2002). Due to 
the peculiar nature of the dispersion relation only waves 
with opposite and equal by magnitude ky can interact, to 
conserve momentum. Thus, only perpendicular cascade 
ensues, which leads to the waves with sufhciently large 
k±, so that nonlinear term is no longer small. It was 
proposed in GS95 that from this point turbulent cascade 
evolves keeping approximate balance between linear and 
nonlinear terms, so called critical balance. The waves 
with larger fey are created by decorrelation in the lat¬ 
eral structure ot the wave packet because of the strong 
interaction. 

In the weak turbulence it is safely to assume that u 
and w wave packets have independent polarization, since 
they interact weakly, only once, and never meet again. 
In strong turbulence this is not nessesarily true as the 
head of the u wavepacket got significatly modified as 
it reaches the tail of the w wave packet. Futhermore, 
as we speak of the Alfven mode, the nonlinear term is 
actually proportional to wk±usin9 where 9 is an angle 
between w and u. Therefore those wave packets that 
have nearly parallel polarizations can survive for longer, 
their cascading being inhibited. 

3. NUMERICAL CODES 

We used several sets of turbulent data flows, produced 
in driven direct MHD numerical simulations. We used 
both incompressible pseudospectral code, which is de¬ 
scribed in detail in GLV02, and compressible isothermal 
ENO code, described in Gho & Lazarian 2002. In both 
cases we used periodic 512^ grid, isotropic random di- 
vergentless driving with correlation time around unity in 
Alfven crossing times. Three incompressible simulations 
had Alfvenic Mach number, calculated as the ratio of the 
kinetic to magnetic energy squared, Ma = 0.7, 1.0 and 
1.4. Gompressible run had Ma = 1.0 and sonic Mach 
number of around unity. 

4. RESULTS 

The sample spectra of the magnetic and kinetic ener¬ 
gies are given in Fig. 1. 

We are asking the question of whether polarizations 
of Alfvenic parts of w and z are truly independent, or 
they have significant correlation. Since the effect we are 
fooking for could be small, we used two different methods 
of separation of scales and modes in order to assure that 
the effect is not due to spurious correlation which came 
e.g. from uncertainty in mode decomposition. 

The first method was to obtain datacubes filtered in 
Fourier space, where only wavevectors around some k±o 
and fcyo were left. The uncertainty in k± was of the order 
of fc_L 0 j while fcyo and its uncertainty was set according 
to GS95 anisotropy. We took the Alfvenic mode which 
has vectors perpendicular to both k and mean B. After 
filtering we tranformed fields into real space. We checked 
that their mean square values were of the order of unity, 
as it has to be for the strong and local turbulence. 

We than proceeded to calculate an rms value of the 


estimated nonlinear term wkj_u \ sind|, the rms of the an¬ 
gular anisotropy, which is < sin^ 0 both compen¬ 

sated by 2^/^, and the estimated nonlinear term without 
sin0, or wk±u. We refer to the ratio of the first and the 
third as the weakening of the interaction corresponding 
to wavevector fc, as this is the number in which nonlinear 
term is smaller than a naive estimate based on indepen¬ 
dency of polarization for w and u. We refer to the second 
term as the measure of geometrical alignment of polar¬ 
ization. It will be unity if w and u are not aligned. 

The second method was one of a transverse structure 
functions (TSF), calculated with respect to the ZocaZ mag¬ 
netic field. Namely, we calculated the square root of the 
forth order structure function 2 < |(fr(;(5usin0|^ >, where 
the S difference is taken between two points, connected 
by vector 1 perpendicular to the local magnetic field, and 
9 is an angle between dw and du, and similar structure 
functions 2 < | sind|^ >, and < \6w6u\‘^ >. And again we 
call the square root of the second TSF the geometrical 
alighnment and the ratio of square roots of the first and 
the third as the weakening of interaction, correspond¬ 
ing to scale 1. Here we replaced full w and z with their 
projections on the plane perpendicular to the B vector. 

The above two methods have different uncertainties 
that go with them. The first method uses global mean 
magnetic field to do both mode decomposition and sepa¬ 
ration of scales. The second method uses local magnetic 
field for separation of scales, as we use 1 perpendicular to 
the local field, but it uses approximate method for mode 
separation, as there is some admix of the slow mode in 
the vectors perpendicular to B. This admix, however, 
has to be small on small scales as the wavevector tends 
to be mostly perpendicular to the magnetic field hence 
the vector of the slow mode is mostly along magnetic 
field. 

If the effects of alighnment are not spurious we will see 
correspondence of two methods when k ~ l/l. We have 
plotted the geometrical alighnment factor and the weak¬ 
ening of the interation factor for some of our datacubes 
on Fig. 2. As we see, there is a correspondence between 
two methods of calculating our factors. 

We see that a purely geometric alighnment factor 
21/2 

< sin^ 9 is quite close to unity. It is certainly 
not enough to produce significant interaction weakening. 
On the other hand, the “weakening factor” is sizable, be¬ 
ing aroung 0.5 for smallest scales. This means that there 
is a significant correlation between wave amplitudes and 
the polarization alignment. In other words, polarization 
alignment is rather strong where there is a large turbu¬ 
lent field. We additionaly confirmed this by plotting the 
distribution of angles between Alfven w and u for re¬ 
gions with different wave amplitudes. It seems that in 
high-amplitude regions the alignment is much stronger. 
Therefore we see that the effect has an intermittent na¬ 
ture. 

5. THE MODEL 

The intermittent weakening of the interaction can in 
principle show itself in different ways. We follow a par¬ 
ticular model that assume, that critical balance is persis¬ 
tent. Indeed, if, in one region, interaction become weak¬ 
ened it is compensated by further growth of k± while fey 
does not grow much because there is no significat lateral 
decorrelation. 



Polarization intermittency 
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Fig. 1.— Spectra of kinetic (solid) and magnetic (dashed) energies of Alfvenic mode. Left: incompressible MHD, = 0.7, hyperdiffusion 
of the third order, lysA^, right: compressible run with Ma = 1, Ms = 1. 
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Fig. 2.— The geometrical alignment (upper curve) and the weakening of interaction factor (lower curve), left panel: incompressible, 
Ma = 0.7, comparison of two methods refered in text, diamonds are for the first method, they were put on distance scale I ~ 0.7/fc, lines 
are for second method (TSF). Solid, dotted and dashed lines are for three simulations, separated by Alfvenic cross time. Central panel: 
incompressible, Ma = 1-0, same notation. Right panel: Compressible, Ma = 1, Ms = 1, we used a global mean magnetic field for both 
mode separation and scale separation. 


Let’s assume that weakening factor scale as Then, 
a critical balance would mean that VA^uk^^ 6u'^k]_ 
Here we use u instead of both u and w, as we consider 
balanced case. Then cascading time could be determined 
by either linear or nonlinear term as r ~ Su~^kL^°‘- 
The Kolmogorov hypothethis du^ jr = const will give us 

the spectrum of Su ~ or Ek ~ 

The anisotropy is obtained from critical balance: fey ~ 

^2/3-2a/3_ 

In case of a = 0, scale independent weakening factor, 
we reproduce GS95 spectrum and anisotropy. In order to 
achieve Iroshnikov-Kraichnan spectrum we have to take 
a = 1/4. However in this case anisotropy will be stronger 
than in GS95, namely fey ~ fc_/ . 

6. INTERMITTENCY 

In section 2 we chose to use 4th order structure func¬ 
tions (SF) to quantify interaction weakening. This choice 
is somewhat arbitrary. Indeed, the Kolmogorov-type ar¬ 
guments suggest using 3rd order functions (see Monin & 
Yaglom 1975), but, alas, we cannot naturally construct 
such a structure function in MHD (for discussion, see 
Biskamp 2003, sec. 7.3.3). If we assume that weakening 
is proportional to the wave amplitude, it is more natural 
to use 4th order SF. We desided to calculate the relative 


scaling exponents for combined Elsasser fields, as well as 
nonlinear term. This will quantify how much depletion 
of interaction or “interaction weakening” that we calcu¬ 
lated, depends on the choice of the order of the SF. 

One of the curious things we have found is that ex¬ 
tended self-similarity (ESS, see Benzi et ah, 1993) be¬ 
tween SEs of different order is much better for nonlinear 
term than for combined elsasser fields. For the latter, 
self-similarity between SFs are mostly limited to iner¬ 
tial range. Obviosly there is no good extended similar¬ 
ity between nonlinear term and elsasser fields. We shall 
cautiosly conclude here, that nonlinear term, raw sin d is 
more “fundamental” than wu. 

We have plotted relative scaling exponents for nonlin¬ 
ear term and the product of Elsasser fields on figure 3. 
We used 4th order SF sin^ 9 as basic, but due to 
ESS it is easy to recalculate them for basic SF of any or¬ 
der. We also plotted a, the weakening factor exponent, 
calculated using SFs of different order. We see that a is 
generally smaller than 1/4. 

7. DISCUSSION 

In our incompressible simulations magnetic field en¬ 
ergy was typically slightly larger that kinetic energy (see 
Fig 1.) This difference is the so-called residual energy 
(see Muller, Grappin, 2005). Residual energy creates a 
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Fig. 3. — Incompressible cubes, Ma = 0.7, Lower panel: 
Relative scaling exponents for SFs < |<5ui(5u|J’/^ > (solid) and 
< |(5ru5usin> (dashed) with p = 4 taken as basic SF. Upper 
panel: interaction weakening scaling factor a (see sec. 5) calculated 
using different order structure functions. 


polarization alignment, but of different type that we mea¬ 
sure. In this case w and u are systematically antiparallel. 
However, we are looking not for difference between p.d.f. 
of angles of 0 and 180 degrees, but for difference for 0 and 
90 degrees. It would be a second-order correction and is 
generally pretty small, as long as kinetic and magnetic 
energies are close. We estimated that in our case this 
spurious effect is never larger than 2%. 

Maron and Goldreich (2001) observed axial asymme¬ 
try or net polarization of Alfven waves. They claimed 
that this effect is strongest in decaying turbulence, but 
in forced turbulence it is bound. They did not, however, 
connected this effect to the flatter than —5/3 spectra 
observed in both decaying and forced runs. They specu¬ 
lated that spectra could be flatter due to intermittency. 

Boldyrev (2005) considers weakening of interaction 
that comes from three-dimensional structure of the eddy. 
As the model he uses same Kolmogorov-type arguments, 
so his formulae can be reproduced by ours by taking 
a = a'/(3 + a'), where a' is Boldyrev’s a. However, 
in our simulations we saw that purely geometrical ar¬ 
guments can not fully describe weakening and the real 
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weakening is due to correlation between angular align¬ 
ment and the wave strength. 

Somewhat unexpected result is that compressible 
cubes show more of a polarization alignment, even 
though they naturally have a smaller inertial range due 
to large numerical diffusivity and viscosity. So far within 
our approach a is not constrained and we do not claim 
that weakening will nessesarily behave as a power-law. 
We leave this two issues to further study. 

While the magnitude of the changes of the spectral in¬ 
dex as well as anisotropy that we observe is not large, 
the potential implications of this can be very substan¬ 
tial. Turbulence in interstellar medium is injected on the 
scale of dozens of parsecs (see Lazarian & Pogosyan 2000, 
Farmer & Goldreich 2004), and on the megaparsec scales 
in clusters of galaxies (see Gassano & Brunetti 2005). 
The scale at which Alfven turbulence dissipates can be 
of the order of thermal proton Larmor radius, i.e. thou¬ 
sands of kilometers. As the result of such a humongous 
scale separation any changes in the spectra and scale- 
dependent anisotropy will have important consequencies 
for the turbulent energy available at a sufficiently small 
scale. Gosmic ray acceleration and propagation are the 
processes that are directly affected. 


8. SUMMARY 

In the paper above we have demonstrated that 

1. Turbulent magnetized flow spontaneously develop 
regions where the mutual shearing of the oppositely mov¬ 
ing Alfven waves is weakened due to polarization align¬ 
ment of the waves. 

2. The amplitude of waves is enhanced in the regions 
of correlated polarization. 

3. Even though the effect of polarization alighnment 
is weak, it affects spectrum and anisotropy, which could 
have a significant impact on a wide range of astrophysical 
phenomena. 
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